function eq = solve_D(param, glob, options, D0)
    % Given a value of C, finds the value of D that
    % makes the free entry condition hold with equality

    D1 = .01*D0;
    D2 = 10*D0;
    
    Dmax = D2;
    
    for t = 1:options.itermax_DC

        param.D = (D1 + D2)/2;
        Ds(t) = param.D;        

        param.w = param.omega*param.C; % HH FOC

        tic
        eq = solve_cL(param,glob,options);     % solve incumbent & entrant pbms 
        eq         = setup_transition_matrices(eq, param,glob,options); % set up transition matrices
        [diff, eq] = find_stationary_dist(eq, param,glob,options);
        toc

        if abs(diff) < options.tolD
            break
        end

        D1s(t) = D1;
        D2s(t) = D2;
        Ds(t) = param.D;

        figure(47)
        subplot(3,1,1)
        hold off
        plot(1:t, D1s, 'rx', 'MarkerSize', 5); drawnow;
        hold on
        plot(1:t, D2s, 'rx', 'MarkerSize', 5); ylim([0 Dmax]); drawnow;
        plot(Ds); title('Guessed value of demand index'); drawnow;

        subplot(3,1,2)
        entry_values(t) = diff;
        plot(1:t, entry_values); title('Difference between D_in and D_out');

        if diff < 0
            D2 = param.D; 
        else
            D1 = param.D;
        end

       
    end
    
    % Also finds the stationary equilibrium
    % First, finding the mass of entrants that 
    % is consistent with the value of D that we found above
        
    figure(8000)
    plot(eq.Lp)
    
end
